{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from collections import defaultdict\n",
    "from math import factorial\n",
    "import sdm as sdmlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "bits = 1000\n",
    "sample = 1000000\n",
    "radius = 451"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "_comb_cache = {}\n",
    "def comb(a, b):\n",
    "    if a < 0:\n",
    "        return 0\n",
    "    if b == 0:\n",
    "        return 1\n",
    "    ret = _comb_cache.get((a, b), None)\n",
    "    if ret is None:\n",
    "        ret = comb(a-1, b) + comb(a-1, b-1)\n",
    "        _comb_cache[(a, b)] = ret\n",
    "    return ret\n",
    "\n",
    "def comb(a, b):\n",
    "    return factorial(a) // factorial(b) // factorial(a-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def genM(n, h):\n",
    "    M = defaultdict(int)\n",
    "    for i in range(n):\n",
    "        x = h+i\n",
    "        y = i\n",
    "        if 0 <= x < n and 0 <= y < n:\n",
    "            M[(x, y)] = comb(n-h, i)\n",
    "\n",
    "        x = h-i\n",
    "        y = i\n",
    "        if 0 <= x < n and 0 <= y < n:\n",
    "            M[(x, y)] = comb(h, i)\n",
    "\n",
    "    for j in range(1, h+1):\n",
    "        mult = M[(h-j, j)]\n",
    "        for i in range(1, n):\n",
    "            x = h+i\n",
    "            y = i\n",
    "            nx = x - j\n",
    "            ny = y + j\n",
    "\n",
    "            if 0 <= x < n and 0 <= y < n and 0 <= nx < n and 0 <= ny < n:\n",
    "                #print('i={} j={} ({}, {}) ({}, {})'.format(i, j, x, y, nx, ny))\n",
    "                M[(nx, ny)] = mult * M[(x, y)]\n",
    "\n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def printM(M):\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            y = n-1-i\n",
    "            x = j\n",
    "            if (M[(x, y)]):\n",
    "                print('{:3d}'.format(M[(x, y)]), end=' ')\n",
    "            else:\n",
    "                print('   ', end=' ')\n",
    "        print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi2_fn(bits, sample, radius, d):\n",
    "    total = 0\n",
    "    M = genM(bits, d)\n",
    "    for i in range(radius+1):\n",
    "        for j in range(radius+1):\n",
    "            total += M[(i, j)]\n",
    "    return 1.0 * sample * total / (2**bits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "_phi_fn_cache = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi_fn(n, H, r, d, steps=500):\n",
    "    key = (n, H, r, d, steps)\n",
    "    if key in _phi_fn_cache:\n",
    "        return _phi_fn_cache[key]\n",
    "    v = []\n",
    "    for _ in range(steps):\n",
    "        bs1 = sdmlib.Bitstring.init_random(n)\n",
    "        bs2 = bs1.copy()\n",
    "        bs2.flip_random_bits(d)\n",
    "        selected1 = address_space.scan_thread2(bs1, r)\n",
    "        selected2 = address_space.scan_thread2(bs2, r)\n",
    "        x = len(set(selected1) & set(selected2))\n",
    "        v.append(x)\n",
    "    mu = 1.0*sum(v)/len(v)\n",
    "    _phi_fn_cache[key] = mu\n",
    "    return mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "address_space = sdmlib.AddressSpace.init_random(bits, sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(280.22, 280.2046153835306, 0.015384616469418688)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = 102\n",
    "a = phi_fn(bits, sample, radius, h, steps=200)\n",
    "b = phi2_fn(bits, sample, radius, h)\n",
    "(a, b, a-b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "72.90908288196275"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1072.25 1071.85004892 0.399951076334\n",
      "964.8 958.436527391 6.36347260943\n",
      "953.9 958.436527391 4.53652739057\n",
      "912.9 902.208104539 10.6918954609\n",
      "901.0 902.208104539 1.20810453914\n",
      "862.8 860.394028924 2.40597107589\n",
      "868.25 860.394028924 7.85597107589\n",
      "826.2 825.8454042 0.354595800216\n",
      "818.25 825.8454042 7.59540419978\n",
      "804.6 795.873633311 8.72636668879\n",
      "788.15 795.873633311 7.72363331121\n",
      "768.2 769.130488339 0.930488339434\n",
      "783.2 769.130488339 14.0695116606\n",
      "742.0 744.827181999 2.82718199914\n",
      "742.35 744.827181999 2.47718199914\n",
      "710.15 722.455125934 12.3051259336\n",
      "722.15 722.455125934 0.305125933626\n",
      "703.95 701.663613943 2.28638605706\n",
      "702.6 701.663613943 0.936386057057\n",
      "678.85 682.198589157 3.34858915679\n",
      "686.85 682.198589157 4.65141084321\n",
      "665.55 663.86892494 1.68107506028\n",
      "671.4 663.86892494 7.53107506028\n",
      "646.8 646.526475663 0.273524336524\n",
      "646.5 646.526475663 0.0264756634762\n",
      "627.65 630.053593212 2.4035932119\n",
      "635.2 630.053593212 5.1464067881\n",
      "614.25 614.354954716 0.104954715998\n",
      "603.55 614.354954716 10.804954716\n",
      "594.35 599.352010301 5.00201030114\n",
      "610.85 599.352010301 11.4979896989\n",
      "589.5 584.979091748 4.52090825168\n",
      "589.5 584.979091748 4.52090825168\n",
      "578.0 571.180612164 6.81938783576\n",
      "571.35 571.180612164 0.169387835762\n",
      "561.5 557.909004388 3.59099561246\n",
      "564.45 557.909004388 6.54099561246\n",
      "548.65 545.123172894 3.52682710641\n",
      "545.9 545.123172894 0.776827106406\n",
      "536.7 532.7873109 3.91268909952\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-1ef1ac5f9df5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mh\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1001\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mphi_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradius\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m     \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mphi2_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradius\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'{} {} {}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-8-ef2a5a8798bb>\u001b[0m in \u001b[0;36mphi_fn\u001b[0;34m(n, H, r, d, steps)\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0mbs2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflip_random_bits\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m         \u001b[0mselected1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0maddress_space\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscan_thread2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbs1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m         \u001b[0mselected2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0maddress_space\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscan_thread2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbs2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mselected1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mselected2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m         \u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/msbrogli/Dropbox/FGV - Doutorado/working-papers/sdm-framework/docs/notebooks/sdm/__init__.pyc\u001b[0m in \u001b[0;36mscan_thread2\u001b[0;34m(self, bs, radius, thread_count)\u001b[0m\n\u001b[1;32m    224\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mscan_thread2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradius\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mthread_count\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    225\u001b[0m         \u001b[0;31m# See https://docs.python.org/3/library/ctypes.html#type-conversions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 226\u001b[0;31m         \u001b[0mselected\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mc_uint\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    227\u001b[0m         \u001b[0mcnt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlibsdm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_scan_thread2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpointer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbs_data\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc_uint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mradius\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mselected\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mthread_count\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    228\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mselected\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mcnt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "for h in range(0, 1001):\n",
    "    a = phi_fn(bits, sample, radius, h, steps=20)\n",
    "    b = phi2_fn(bits, sample, radius, h)\n",
    "    print('{} {} {}'.format(a, b, abs(a-b)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
